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INTRODUCTION 


Grid smoothing and orthogonal ization procedures have been developed 
and implemented in the construction of two and three-dimensional grids. 
The procedures are based on the variational methods of grid generation. 
The results of this work are included in the attached paper which will 
be presented at the AIAA Aerospace Sciences Meeting in Reno. The 
two-dimensional examples were computed using the MSU IRIS Graphics 
Workstation. Color plots, which could not be included in the paper, 
will be presented at the meeting. The attached paper also contains some 
three-dimensional examples which are the product of some work done 
during the summer at the Air Force's Arnold Engineering Development 
Center. That work gave us the opportunity to test our methods on large 
three-dimensional grids without using any funds from the Grant. 

There is one result of our work during this reporting period which 
is most significant. It has been demonstrated that the elliptic grid 
generation equations, with arbitrary forcing functions, can be solved, 
in their variational formulation, using a gradient method. Since 
gradient methods have a global convergence property, the divergence 
problems often encountered when using SOR iterative methods can be 
avoided. It is not to be concluded, however, that SOR methods should be 
abandoned, since gradient methods tend to converge very slowly. In 
fact, slow convergence was the major problem encountered in the 
three-dimensional grids. 

Further progress has been made on the continuing effort to develop 
conservative interpolation formulas for overlapping grids. The basic 



algorithm was derived earlier and included in our previous semiannual 
report. Since then, our work has been directed toward solving some of 
the technical problems such as: how to efficiently find grid 

intersection points, how to store interpolation coefficients, and how to 
implement the extrapolation procedures at outflow points on the overlap 
boundary. Algorithms have been developed for solving these problems. 
Final verification of the algorithms has not been completed at this 
time. This work will continue into the next contract period and final 
results will be included in a future report. 
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Abstract 

Optimization methods are developed and ana- 
lyzed for the solution of all the popular varia- 
tional problems in grid generation. Comparisons 
are made between the different variational methods. 


I. Introduction 

Variational methods in grid generation have 
recently become a viable alternative’ to the more 
established elliptic methods. Their appeal is 
likely due to the geometric foundation of the 
variational methods. These methods were derived 
using grid properties such as distances between 
grid points, orthogonality of grid lines, and the 
volumes of grid cells. Because of the earlier 
successes of elliptic methods, the first solution 
algorithms consisted of converting the variational 
problem to a problem of solving the elliptic Euler 
equations. The Euler equations were then solved 
using available finite difference algorithms. This 
approach was used by Brackbill and Saltzman [1], 
and Roache and Steinberg {4]. One disadvantage in 
this procedure is that the Euler equations are 
considerably more complicated than the original 
variational integral. Due to the complexity of the 
Euler equations, they are difficult to solve and 
solution algorithms may not converge. These diffi- 
culties have led to the development of algorithms 
for the direct solution of the variational problem. 
This approach was followed by Kennon and 
Dulikravich [3] and Carcaillet et al. [2]. In 
these two reports, different discrete problems were 
formulated as unconstrained optimization problems 
and then solved using a conjugate gradient itera- 
tive method. There are several advantages in this 
alternative. The original geometric properties of 
the grid are stated precisely by the discrete 
problem. Thus there is no truncation error as 
there is in the numerical solution of the Euler 
equations. Also, the conjugate gradient and simi- 
lar descent methods will always converge to a local 
extrema. However, the convergence rate may be very 
slow. 
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This report will present a survey of the types 
of integrals that may be included in a variational 
problem and the geometric properties that each 
integral imposes upon the grid. Optimization 
methods for the direct solution of the variational 
problems will also be covered. A Jacobi-Newton 
iterative method will be developed and compared 
with the Fletcher-Reeves conjugate gradient method. 
By combining the simpler Jacobi-Newton iteration 
with the direct solution of the variational prob- 
lem, one has an algorithm which is easier to pro- 
gram than the Euler equations and requires less 
storage and less computer time per iteration than 
the conjugate gradient method. 


II. Variational Integrals 

Variational methods may be used to influence 
three geometric properties of a grid. These are 
the smoothness of the grid point distribution, the 
size of the grid cells, and the orthogonality of 
the grid lines. Each property, smoothness, cell 
volumes, and orthogonality, can be assigned a 
numerical value by evaluating an integral. Most of 
the results will deal with three-dimensional prob- 
lems, and the physical and computational variables 
will be denoted by x, y, z and £ , n, C. respec- 
tively. 

There are two basic smoothness integrals, one 
using the gradient of the physical variables and 
the other using the gradient of the computational 
variables. The integrals are denoted as and 

with D and R the computational and physical regions. 

S 1 = j D p2!r J 2 + q2 ! r n 1 2 + r2 I | 2 dedndC (1) 

S- *= ( P 2 |vc | 2 + Q 2 1 Vn 1 2 + R 2 1 V£ | 2 dxdydz 
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where J is the Jacobian and r » (x,v,z). 
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The functions P, Q, and R, may be selected to con- 
trol the grid point distribution. They directly 
affect the grid aspect ratio, and their influence 
can be readily seen by considering a one- 
dimensional transformation of any of the computa- 
tional variables. The functions minimizing either 
of these integrals will define a smooth mapping 
between the physical and computational regions in 
the sense that the mapping functions are harmonic 
and therefore infinitely differentiable. However, 
the mapping defined by the minimization of may 

not be one-to-one. The mapping defined by the 
minimization of is one-to-one, but the con- 
struction of the mapping is more difficult. 
Current methods require the numerical solution of 
a system of nonlinear Euler equations. 

For the two-dimensional case, there is a sim- 
ple relation between the above integrals, (1) and 
(2), and the elliptic equations commonly used for 
grid generation. A two-dimensional version of 
in (2) may be written in the form 


r q 2 |vc| 2 + p 2 |v n | 2 dxdy 


- f P 2 |r,| 2 + Q 2 (r | 2 )/J dWn. 
J D n 


where J is the Jacobian 


The Euler 


3(5, n)' 

equations for minimizing this integral may be 
written as 


A uniform distribution of grid cell volumes 
vould not be appropriate in most practical appli- 
cations. However, a near constant value of some 
weighted volume is often called for. The weight 
may depend on some desired distribution of volumes 
or a solution-dependent quantity used to define an 
adaptive grid. From a variational point of view, 
the optimal volume distribution could be defined 
by minimizing the integral 


V 
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Gdxdydz 


(3) 

j GJd£dndC 
J D 


where G » Jw, w ■ w(5,n,C)- 

While there is no guarantee that the variational 
problem will have a solution which gives a non- 
singular mapping between physical and computa- 
tional spaces, it can be noted that if such a 
solution exists, then the weighted Jacobian G will 
be constant. This follows from the Euler 
equations. The argument is simpler in the two- 
dimensional case. The Euler equations are then 


Vn * Vs 
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- G r x + G x r * 0, 

S n n S 
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and if it is assumed that the Jacobian 
x^y^ - x^y^ #0, it follows that 


v 2 n « -UnPyvni 2 


G - G - 0 
S n 


or in computational space, with r * as 


or G must be constant. The constant is determined 
by the physical region and the value of w since 
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There is one other volume integral that should 
be considered. If 


This is the elliptic system with control functions 


4 m ( f nQ ) ^ , ip - UnP)^. 

If # and are given, then the corresponding 
control functions for the variational integral are 


P « exp( Udn)*Q * exp( UdS). 


V - f |VG| 2 dedn<iC, (5) 

J D 

then the integrand is a second order differential 
expression in the physical variables x, y, and z. 
Consequently, the Euler equations are a system of 
fourth order equations. Now there must be two 
boundary conditions if a unique solution to the 
Euler equations is to be determined. One condition 
is imposed by fixing the boundary correspondence. 
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A second condition can be motivated by considering 
a lower order problem where G, the integrand of 
, is a single dependent function of the variables 

Ht and C- The Euler equation would then be 
Laplaces equation 


( p-"- 

J D \lM IK, 





v 2 g - 0. 


The equation has a unique solution if G is given 
on the boundary of the computational region. The 
function G would be uniquely determined up to a 
constant value if on the boundary. 


G^ - 0 for £ * constant 
G - 0 for n ■ constant 

n 

G^ « 0 for C * constant. 



Irrespective of the solution algorithm selected, 
the minimization of A ^ is a more difficult problem 

since the integrand is a rational function of the 
physical derivatives rather than a multinomial 
as in A^. 

With the variational grid generation methods 
one selects a nonnegative linear combination of the 
form 


The constant value would again be determined by 
(4). Either of these types of boundary conditions 
could be advantageous in grid generation. Being 
able to specify the value of G along the boundary 
would aid in maintaining desired distributions of 
grid points near boundaries such as is needed in 
grids for the solution of problems with boundary 
layers. On the other hand, the derivative bound- 
ary conditions for G would be of value when 
several grids are patched together and it is 
desired that the cell sizes vary continuously from 
one subregion into another. The latter condition 
could also be used to generate a more uniform dis- 
tribution of G over a single region. 

If a grid is orthogonal, then 


I * cxS i + + yA k , (8) 

where i, j, and k are either 1 or 2. Since the 

S . , V., and A. may be of differing magnitudes, the 
l 3 K 

coefficients a, 3, and y may contain scaling 
factors. One method of scaling that can be used is 

to compute values V^ 0 ^, and A^^ from an 

initial algebraic grid and then select non-negative 
numbers a, b, and c so that 

a + b + c * 1. 
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where r * (x,y,z). Therefore, one measure of 
skewness is the integral 


Now if we choose 
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) 2 + (r *r ) 2 + (r *r J 2 d£dndC. (6) 
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then the initial value of I is the same for any 
combination of a, b, and c. 


The integrand depends on the size of the grid 

cells since, for example. III. OPTIMIZATION 


r e • r n ■ l l r e l l I l r n l I cos 6 

where 8 is the angle of intersection between the 
U ana n grid lines. As a result, minimizing the 
integral A^ would have little effect in subregions 

with extremely small cells. For grids with widely 
varying cell volumes, the following integral would 
be more appropriate 


In the direct solution of the variational 
problem each integral is approximated on a square 
grid in computational space. Thus the continuous 
problem of minimizing I is reduced to the discrete 
problem of finding a minimum value of an objective 
function 

F ‘K J.. <»> 
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where i, j f k ranges over the indices of the grid 

points. The grid points on the boundary would be 

included although it is assumed that these points 

are fixed. Each term F. . , is to contain all 
i • j > k 

difference approximations centered at the grid 
point r which are used in the formation of 

1 » j 

the objective function F. The optimal grid is 
found by obtaining the unconstrained minimizer of 
the objection function F. There are many algo- 
rithms for solving optimization problems. An 
algorithm which has been successfully implemented 
in grid generation is the Fletcher-Reeves conju- 
gate gradient method. The method can be used with 
any of the variational integrals of the previous 
section. Since the approximations of (2) and (7) 
involve quotients, it is possible that somewhere 
in the iteration scheme a division by a number 
close to zero may occur. This problem can be 
easily avoided by placing a lover bound on the 
denominators. Only a two-dimensional grid prob- 
lem using (2) has been attempted because the three- 
dimensional integral becomes very complicated when 
transformed to computational space.. 

The conjugate gradient and related descent 
methods have the attractive feature of always con- 
verging to a grid which represents a local minimum 
value for the objective function. Their use in 
the solution of large three-dimensional grid 
generation problems may not be very efficient. At 
each iteration a one-dimensional minimization 
problem must be solved to determine the optimal 
step-size in the descent direction. The step-size 
problem can not be easily solved and frequently 
requires several evaluations of the objective 
function. This is especially true when the inte- 
grals in (2) and (7) are included. Therefore, 
alternate solution algorithms have been investi- 
gated. The application of Newton or quasi-Newton 
methods for nonlinear minimization problems in 
grid generation does not appear to be practical 
because of the number of matrix operations. 
Motivated by algorithms which are used to solve 
the discretized Euler equations, a Jacobi-Newton 
method has been selected as an alternate algorithm 
for solving the optimization problem. The 
following development covers the basic steps of 
the algorithm. 


The system of nonlinear equations is solved using 
a Jacobi-Newton iterative method. At each point, 
Newton's method is used to approximately solve the 
three equations given in (10). In all of our 
examples, only one Newton iteration was performed. 
In fact, for certain variational integrals, the 
system is linear and the exact solution is obtained 
in one iteration. 

The Jacobi-Newton iterative method described 
above has been used to solve a wide range of two 
and three-dimensional grid generation problems. 

So far, it has not been used with the variational 
integrals and A^ because the quotients would 

cause the Newton iteration equations to be very 
lengthy. Convergence problems have not been 
encountered except in some cases where the grid 
folded. It is often necessary to use under- 
relaxation for larger problems. 


IV. RESULTS 


The following examples are presented to 
demonstrate the effect of the grid optimization 
procedure when applied to a large-scale three- 
dimensional grid. As can be seen, the success of 
the optimization process depends primarily on the 
proper choice of the variational integrals included 
in (8). When a decision had to be made between 
efficiency and robustness, it was the policy to 
develop a robust algorithm. A numerical step-size 
determination was used in the conjugate gradient 
iteration even though it was sometimes possible to 
calculate the exact optimal step-size. 

All of the variational integrals have been 
used in the optimization of small two-dimensional 
grids. In all cases the objective was to improve 
an algebraic grid while keeping the basic grid 
point distribution. Thus the functions P and Q of 
integrals and S 2 are defined simply as 


P 
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l r (o)l 
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, Q 
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A local minimum value of F occurs at points 

where 


and the weight function w 


in and 


is 


VF - 0, 

or writing this system in component form 


w(e.n) - ( |r (0) | 
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dx 
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where the superscript (0) is used to indicate that 
these functions are computed from tne coordinates 
of the initial algebraic grid. One-sided differ- 
ences have been used in the approximation of all 
variational integrals. Consequently, at each grid 
point there are four possible approximations for 
each integrand. In order to retain any grid 
symmetry, the average of these four values is used 
as the approximation of the integrand. 

The integral generates the only optimal grid 
which is smooth and unfolded in the neighborhood 
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of concave corners. Such a grid is plotted in 
Figure 1. The volume integrals and would 

not generate smooth grids and the integral 

would generate a grid which folds over the corner. 
However* does generate acceptable grids near 

convex corners. For roost regions, the orthogo- 
nality integrals, and A 2> cannot be used alone 

to solve grid generation problems. An improved 
grid can often be generated by choosing the 
objective function as a convex combination of an 
orthogonality integral and either a smoothness or 
volume integral. The second integral is used to 
prevent grid folding. Figures 2 and 3 demonstrate 
the characteristics of the two orthogonality 
integrals. The grid in Figure 2 is generated 
using A^ and S^. Note the effect of including 

A^. , While there is some decrease in skewness 

where the grid is coarse, there are other regions 
where the grid lines are even more skewed than in 
the original interpolated grid. That is not the 
case when A£ and are used to form the objective 

function. That grid is plotted in Figure 3. In 
order to emphasize the effect of the orthogonality 
integrals in the last two examples, only enough 
smoothness was included to prevent grid lines from 
crossing. The volume integral has been very 

helpful in the generation of grids near singular- 
ities. An example of a grid generated by minimiz- 
ing V is presented in Figure 4. Figures 5 and 6 
give two more examples where the grid optimization 
process can reduce skewness in a grid without 
destroying the grid point distribution. The 
integral was generally less effective in these 

examples. Convergence was also slower, as might 
be expected with a method involving higher order 
derivatives. 

Most of the grid optimization options have 
been incorporated in a state-of-the-art grid 
generation code developed at AEDC. The following 
examples illustrate the optimization of large 
three-dimensional grids for practical geometric 
configurations. Figure 7 represents the grid for 
an internal flow problem. Figure 8 gives a view 
of a part of a grid surface after optimization. 
Figure 9 is the grid for an external flow problem. 
The effect of the grid optimization process is 
illustrated in Figure 10. 



Figure 1. Optimized grid using $ 2 


It should be mentioned that the term "volume'’ 
used in this report is actually a misnomer. The 
Jacobian in equations (3) and (4) should always be 
computed from the definition and the appropriate 
difference approximations. If the cell areas (or 
volumes) are substituted for the Jacobians, then 
the optimal grid may contain triangular or even 
non-convex cells. 


V. CONCLUSIONS 

The successful application of variational 
methods in grid optimization depends on tailoring 
the integral to be minimized to the shape of the 
physical region and the desired distribution of 
grid points. If an algebraically generated grid 
is to be optimized, control functions can be 
defined using the arc length distribution of points 
along each family of grid lines. The solution of 
the variational problem can be computed using the 
traditional optimization methods or by iterative 
methods similar to those currently being used to 
solve elliptic partial differential equations. 

The convergence rate for either iterative method is 
generally slow when applied to large three-dimen- 
sional grids. 
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Figure 2. Optimized grid using S-j and A-j 
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